rm(list = ls(all = TRUE))
setwd("C:\\Users\\lenovo\\Desktop\\R")

#主成分分析（Principal Component Analysis， PCA）
#主成分分析是在确保数据信息损失最小的原则下，
#把多个指标转化为是少数几个不相关的综合指标的数据降维关系

#因子分析（Factor Analysis，FA）
#因子分析通过寻找一组更小的、潜在的或者隐藏的结构
#来解释已经观测到的变量之间的关系

#12.1 主成分分析
#Hotelling，1933，本质为降维
#12.1.1 主成分的定义

#本质就是对原始变量进行线性变换
#主成分就是原始变量的线性组合
#各线性组合的组合系数均为单位向量
#在实际应用中，所确定的主成分的个数总是小于原始指标的个数
#实际只有数据中存在多个相关变量时才需要用到主成分分析

#12.2.2 主成分的求解

#本质是对自变量的协方差矩阵（或相关系数矩阵）进行特征分解
#即，求其特征值和特征向量
#特征值代表了特征的重要程度，特征向量代表了坐标旋转后的特征方向
#base包中的eigen

library(datasets)
cor.matrix <- Harman23.cor$cov
eigen(cor.matrix)

#主成分的个数遵循的两大原则：
#前k个主成分的累计贡献率达到某一特定的值（如70%-80%）
#保留的主成分对应的特征值大于1（Kaiser-Harris准则）
#实际问题中还需要结合问题的背景和主成分的实际意义

eigen(cor.matrix)$values
#结果表明，前两个特征值大于1，
#且前两个特征值之和占总的特征值之和的80.5%
#接下来由特征值对应的单位向量获得各个主成分（查看时保留三位小数）
round(eigen(cor.matrix)$vectors,3)

#[,1]和[,2]是自变量的系数
#第一个主成分与每个指标都为负相关，可以看作一个综合指标
#第二个主成分，前4个指标（长度）负相关，后4个指标（体型）负相关，
#可以看作一个体型指标（结合实际意义）

#上述计算利用stats包里的函数princomp轻松完成
#可以使用原始数据或者原始数据的协方差矩阵（或相关系数矩阵）
#本例为相关系数矩阵
PCA <- princomp(covmat = cor.matrix)
summary(PCA, loadings = TRUE)
#对于summary，参数loadings默认为False，设为True显示载荷矩阵
#载荷矩阵中，绝对值小于0.1的数据不显示，
#若要显示，参数cutoff设为0

#借助点线图辅助判断
#碎石图
#scree一词来自地质学，表示在岩层斜坡下方发现的小碎石
#这些小碎石地质学价值不高，可以忽略
screeplot(PCA, type = "lines")
abline(h = 1)
#前2个主成分的散点在斜坡上
#后6个主成分的散点形成了平台，且特征值均小于1
#因此保留前2个主成分

#为了分析主成分的意义，提取载荷矩阵
load <- loadings(PCA)
plot(load[,1:2], xlim = c(-0.5,0.5), ylim = c(-0.5,0.5))
#图中点打上标签
text(load[,1],load[,2],adj = c(-0.3,0))
#加上线
abline(h = 0, v = 0)
#进一步表明了两个主成分具有的前面已经讨论过的特征

#12.2.3 主成分分析的注意事项
#主成分分析不会增加总信息量
#当原始变量之间的相关性较小时，应用主成分分析是没有意义的
#主成分分析可以基于协方差矩阵或者相关系数矩阵，但结果有所不同
#指标之间的取值范围彼此相差不大时可以用协方差矩阵，
#尽量保留原始变量的实际意义
#指标间取值范围相差较大或者量纲不同是，采用相关系数矩阵

#主成分分析的目的是降维，很多时候并不是最终目的

#12.2 因子分析
#20世纪，Charles Spearman和Karl Pearson
#一种主成分分析的推广，
#根据相关性把变量分组，使得组内变量的相关性较高
#但不同组的变量之间的相关性较低
#则每组变量可代表一个基本结构，称为因子
#因子反映已经观测到的相关性
#各公共因子的均值为0，方差为1，且相互独立

#12.2.1 因子分析模型的定义
#指标=公共因子*载荷+特殊因子

#12.2.2 因子分析模型的求解
#原始变量相关矩阵-特殊因子相关矩阵=约化相关矩阵
#约化相关矩阵=载荷矩阵*载荷矩阵的转置矩阵
#其对角线上的元素在0-1之间，称为公共度
#公共度为1-特殊因子的方差
#接近于0，则原始指标受公共因子影响不大
#等于1，则原始指标可由公共因子的线性组合表示，与特殊因子无关
#为了求解载荷，需要对公共因子有所估计，
#估计方法不同，所用的因子分析的方法也就不同
#常用的估计约化相关矩阵的方法为主成分法、主因子法、极大似然法

cormat <- matrix(c(1,0.439,0.410,0.288,0.329,0.248,
                   0.439,1,0.354,0.354,0.320,0.329,
                   0.410,0.354,1,0.164,0.190,0.181,
                   0.288,0.354,0.164,1,0.595,0.470,
                   0.329,0.320,0.190,0.595,1,0.464,
                   0.248,0.329,0.181,0.470,0.464,1),
                 nrow = 6)

#虽然stats包中的factanal可用于因子分析，
#但推荐使用psych包的一系列函数
library(psych)
#第一个参数：原始数据库或者相关系数矩阵
#第二个参数：当给出相关系数矩阵时，提供观测的个数
#第三个参数：指定提取公共因子的方法
#默认为最小残差法（minres），现使用极大似然法（ml）
fa.parallel(cormat, n.obs = 220, fm = "ml")

#结果表明，可以选择保留一个主成分（碎石检验、平行检验）
#或者两个（特征值大于1）

#如果使用因子分析，需要提取两个因子
#对于因子分析，Kaiser-Harris原则要求选择特征值大于0的因子

#使用fa提取公共因子并求解因子载荷矩阵
#提取因子有最小残差法（minres）、主成分法（pa）、极大似然法（ml）
#加权最小二乘法（wls）、广义加权最小二乘法（gls）
FA1 <- fa(cormat, nfactors = 2, n.obs = 220, rotate = "none", fm = "ml")
FA1

#结果显示，第2、4、5、6这4个变量在第一个因子上的载荷较大
#分别为0.57、0.74、0.72、0.60
#另外两个变量（1、3）在两个因子上的载荷比较平均
#因子载荷的意义较难解释，此时使用因子旋转有助于解释

#12.2.3 因子旋转

#因子分析的主要目的不仅仅是找出公共因子，
#更重要的是要弄清公共因子的专业意义
#数学上可证明，满足因子分析模型的公共因子不唯一
#只要旋转初始公共因子，就可以获得一组新的公共因子

#旋转的方法有正交旋转（orthogonal rotation）和斜交旋转（oblique rotation）
#正交旋转：保持各指标的共性方差不变，且旋转后各因子之间仍然互不相关
#斜交旋转：不保证旋转后的各公共因子互不相关，且对因子载荷的解释要复杂得多
#但斜交旋转取得的效果一般优于正交旋转
#在实际分析中，可以先用正交旋转，不易解释时，再用斜交旋转

#在正交旋转中，最常用的是方差最大旋转（varimax）
#其他方法还有四次方最大旋转（quartimax），均方最大选择（eequmax）

#用方差最大旋转重新建立模型
FA2 <- fa(cormat, nfactors = 2, n.obs = 220, rotate = "varimax", fm = "ml")
FA2

#结果显示，公共因子更易解释了
#第一个因子的后3个变量（算术、代数、几何）载荷较大
#可称第一个因子为数学因子
#第二个因子的前3个变量（盖尔语、英语、历史）载荷较大
#可称第二个因子为文史因子

#利用函数fa.diagram生成图形直观地展示因子分析的结果
#digits为小数位数，默认为1位
fa.diagram(FA2,digits = 2)

#正交旋转强制两个因子不相关，如果想允许两个因子相关，可以使用斜交旋转
#在斜交旋转中，最常用的两种方法是oblimin旋转和promax旋转
#install.packages("GPArotation")
FA3 <- fa(cormat, nfactors = 2, n.obs = 220, rotate = "promax", fm = "ml")
FA3

# With factor correlations of 
#ML1  ML2
#ML1 1.00 0.56
#ML2 0.56 1.00

#旋转后因子结构更为明确，同样可以得到一个数学因子与文史因子
#因子相关矩阵显示两个因子的相关系数为0.56，相关性较强
#如果因子之间的相关性很低，可能需要使用正交旋转来重新简化问题
fa.diagram(FA3, digits = 2)

#12.2.4 因子分析的注意事项
#1 因子分析的解不唯一
#综合考虑；结合专业意义
#2 因子得分问题
#相对于主成分分析，因子分析并不是很关注因子得分的计算，
#如果使用原始数据，在函数fa中设置参数score为TRUE即可获得因子得分
#主成分分析的得分是精确的
#但是因子分析的得分是估计出来的
#估计方法有回归法、最小二乘法等

#3 主成分分析与因子分析的关系
#主成分是一种线性变换，对于每一个原始数据矩阵框而言，
#其主成分系数矩阵是唯一的，
#因子旋转用于对因子载荷进一步简化，使得各公共因子具有明确的实际意义
#当特殊因子的变差为0时，因子分析与主成分分析完全等价

#主成分分析侧重于综合原始变量的信息，
#因子分析侧重于解释原始变量之间的关系

#4 探索性因子分析（EFA）与验证性因子分析（CFA）
#探索性：不知道因子个数
#验证性：研究者根据专业知识与经验对潜在因子的结构已有认知，
#提出观测变量与潜在因子之间存在某种假设的关系
#通过现有资料来验证这种假设，并评价观测变量和因子之间的联系
#EFA: Exploratory Factor Analysis
#CFA: Confirmatory Factor Analysis
#CFA属于结构方程模型（Structural Equation Model，SEM）的一种
#install.packages("lavaan")
#install.packages("sem")

??sem
